module gauss_leg_mod
    implicit none
    public:: gauss_leg,func1d
    private
    real(kind=8),parameter::xi5(5)=[-0.90617984593866399281d0, -0.53846931010568309105d0,&
        & 0.0d0, 0.53846931010568309105d0, 0.90617984593866399281d0]
    real(kind=8),parameter::wi5(5)=[0.23692688505618908751d0, 0.47862867049936646804d0,&
        & 0.56888888888888888888d0, 0.47862867049936646804d0, 0.23692688505618908751d0]
    real(kind=8),parameter::xi6(6)=[-0.93246951420315202781d0, -0.66120938646626451366d0,&
        & -0.23861918608319690863d0, 0.23861918608319690863d0, 0.66120938646626451366d0, 0.93246951420315202781d0]
    real(kind=8),parameter::wi6(6)=[0.17132449237917034504d0, 0.36076157304813860757d0,&
        & 0.46791393457269104739d0, 0.46791393457269104739d0, 0.36076157304813860757d0, 0.17132449237917034504d0]
    abstract interface
        function func1d(x) result(res)
            real(kind=8),intent(in)::x
            real(kind=8)::res
        end function func1d
    end interface
    real(8),parameter::eps=1.0d-14
contains
    recursive function gauss_leg(a,b,f)result(res)
        real(kind=8),intent(in)::a,b
        procedure(func1d) :: f
        real(kind=8)::res
        real(kind=8)::xm
        real(kind=8)::res5,res6
        xm=0.5d0*(a+b)
        res5=gauss_leg_node5(a,xm,f)
        res6=gauss_leg_node6(a,xm,f)
        if(abs(res5-res6)<eps)then
            res=res6
        else
            res=gauss_leg(a,xm,f)
        end if
        res5=gauss_leg_node5(xm,b,f)
        res6=gauss_leg_node6(xm,b,f)
        if(abs(res5-res6)<eps)then
            res=res+res6
        else
            res=res+gauss_leg(xm,b,f)
        end if
    end function gauss_leg

    function gauss_leg_node5(a,b,f)result(res)
        real(kind=8),intent(in)::a,b
        procedure(func1d) :: f
        real(kind=8)::xi,xm,res,w
        integer::i
        xi=0.5d0*(a+b)
        xm=0.5d0*(b-a)
        res=0.0d0
        do i=1,5
            w=wi5(i)*f(xi+xi5(i)*xm)
            res=res+w
        end do
        res=res*xm
    end function gauss_leg_node5

    function gauss_leg_node6(a,b,f)result(res)
        real(kind=8),intent(in)::a,b
        procedure(func1d) :: f
        real(kind=8)::xi,xm,res,w
        integer::i
        xi=0.5d0*(a+b)
        xm=0.5d0*(b-a)
        res=0.0d0
        do i=1,6
            w=wi6(i)*f(xi+xi6(i)*xm)
            res=res+w
        end do
        res=res*xm
    end function gauss_leg_node6
end module gauss_leg_mod

program main
    use gauss_leg_mod
    implicit none
    real(kind=8)::res
    real(kind=8)::a,b
    a=0.d0
    b=10.d0
    res=gauss_leg(a,b,func)
    write(*,*)"num  ",res
    write(*,*)"exact",0.2d0*atan(2.d0)
contains
    function func(x) result(res)
        real(kind=8),intent(in)::x
        real(kind=8)::res
        res=1.0d0/(x*x+25.d0)
    end function func
end program main
